{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2.8 Discontinuous Galerkin Methods\n",
    "\n",
    "* Use discontinuous finite element spaces to solve PDEs. \n",
    "* Allows upwind-stabilization for convection-dominated problems\n",
    "* Requires additional jump terms for consistency \n",
    "\n",
    "Interior penalty DG form for $-\\Delta u$:\n",
    "\n",
    "$$\n",
    "\\DeclareMathOperator{\\Div}{div}\n",
    "A(u,v) = \\sum_T \\int_T \\nabla u \\nabla v\n",
    "-  \\sum_F \\int_F \\{ n \\nabla u \\} [v] \n",
    "-  \\sum_F \\int_F \\{ n \\nabla v \\} [u] \n",
    "+ \\frac{\\alpha p^2}{h} \\sum_F \\int_F [u][v]\n",
    "$$\n",
    "\n",
    "with jump-term over facets:\n",
    "$$\n",
    "[u] = u_{left} - u_{right}\n",
    "$$\n",
    "\n",
    "and averaging operator\n",
    "$$\n",
    "\\{ n \\nabla u \\} = \\tfrac{1}{2} (n_{left} \\nabla u_{left} + n_{left} \\nabla u_{right})\n",
    "$$\n",
    "\n",
    "DG form for $\\Div (b u)$, where $b$ is the given wind:\n",
    "\n",
    "$$\n",
    "B(u,v) = -\\sum_T b u \\nabla v + \\sum_F \\int_F b\\cdot n   u^{upwind} v \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import netgen.gui\n",
    "%gui tk\n",
    "from netgen.geom2d import unit_square\n",
    "from ngsolve import *\n",
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The space is responsible for allocating the matrix graph. Tell it that it should reserve entries for the coupling terms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "order=4\n",
    "fes = L2(mesh, order=order, dgjumps=True)\n",
    "u,v = fes.TnT()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Every facet has a master element. The value from the other element is referred to via the\n",
    "`Other()` operator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "jump_u = u-u.Other()\n",
    "jump_v = v-v.Other()\n",
    "n = specialcf.normal(2)\n",
    "mean_dudn = 0.5*n * (grad(u)+grad(u.Other()))\n",
    "mean_dvdn = 0.5*n * (grad(v)+grad(v.Other()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Integrals on facets are computed by setting `skeleton=True`. This iterates over all internal facets. Additionally setting `BND` iterates only over boundary facets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 4\n",
    "h = specialcf.mesh_size\n",
    "a = BilinearForm(fes)\n",
    "a += SymbolicBFI(grad(u)*grad(v))\n",
    "a += SymbolicBFI(alpha*order**2/h*jump_u*jump_v, skeleton=True)\n",
    "a += SymbolicBFI(alpha*order**2/h*u*v, BND, skeleton=True)\n",
    "a += SymbolicBFI(-mean_dudn*jump_v -mean_dvdn*jump_u, skeleton=True)\n",
    "a += SymbolicBFI(-n*grad(u)*v-n*grad(v)*u, BND, skeleton=True)\n",
    "a.Assemble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = LinearForm(fes)\n",
    "f += SymbolicLFI(1*v)\n",
    "f.Assemble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes, name=\"uDG\")\n",
    "gfu.vec.data = a.mat.Inverse() * f.vec\n",
    "Draw (gfu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "DG requires a lot of additional matrix entries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a nze: 23400\n",
      "a2 nze: 6750\n"
     ]
    }
   ],
   "source": [
    "print (\"a nze:\", a.mat.nze)\n",
    "fes2 = L2(mesh, order=order)\n",
    "a2 = BilinearForm(fes2)\n",
    "a2 += SymbolicBFI(u*v)\n",
    "a2.Assemble()\n",
    "print (\"a2 nze:\", a2.mat.nze)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we are solving a convection-diffusion problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 4\n",
    "h = specialcf.mesh_size\n",
    "acd = BilinearForm(fes)\n",
    "acd += SymbolicBFI(grad(u)*grad(v))\n",
    "acd += SymbolicBFI(alpha*order**2/h*jump_u*jump_v, skeleton=True)\n",
    "acd += SymbolicBFI(alpha*order**2/h*u*v, BND, skeleton=True)\n",
    "acd += SymbolicBFI(-mean_dudn*jump_v -mean_dvdn*jump_u, skeleton=True)\n",
    "acd += SymbolicBFI(-n*grad(u)*v-n*grad(v)*u, BND, skeleton=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `IfPos` checks whether the first argument is positive. Then it returns the second one, else the third one. This is used to define the upwind flux. The check is performed in every integration-point on the skeleton:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = CoefficientFunction( (20,1) )\n",
    "acd += SymbolicBFI(-b * u * grad(v))\n",
    "uup = IfPos(b*n, u, u.Other())\n",
    "acd += SymbolicBFI(b*n*uup*jump_v, skeleton=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "acd.Assemble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes)\n",
    "gfu.vec.data = acd.mat.Inverse(freedofs=fes.FreeDofs(),inverse=\"umfpack\") * f.vec\n",
    "Draw (gfu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hybrid Discontinuous Galerkin methods\n",
    "use additionally the *hybrid* facet variable on the skeleton:\n",
    "\n",
    "$$\n",
    "\\DeclareMathOperator{\\Div}{div}\n",
    "A(u,\\widehat u; v, \\widehat v) = \n",
    "  \\sum_T \\int_T \\nabla u \\nabla v\n",
    "- \\sum_T \\int_{\\partial T} n \\nabla u (v-\\widehat v)\n",
    "- \\sum_T \\int_{\\partial T} n \\nabla u (u-\\widehat u)\n",
    "+ \\frac{\\alpha p^2}{h} \\sum_F \\int_F (u-\\widehat u)(v-\\widehat v)\n",
    "$$\n",
    "\n",
    "the jump-term is now replaced by the difference $u - \\widehat u$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No additional matrix entries across elements are produced. Dirichlet boundary conditions are set as usual to the facet variable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "order=4\n",
    "V = L2(mesh, order=order)\n",
    "F = FacetFESpace(mesh, order=order, dirichlet=\"bottom|left|right|top\")\n",
    "fes = FESpace([V,F])\n",
    "u,uhat = fes.TrialFunction()\n",
    "v,vhat = fes.TestFunction()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, the jump is the difference between element-term and facet-term:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "jump_u = u-uhat\n",
    "jump_v = v-vhat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A non-zero elements: 5825\n"
     ]
    }
   ],
   "source": [
    "alpha = 4\n",
    "condense = True\n",
    "h = specialcf.mesh_size\n",
    "n = specialcf.normal(mesh.dim)\n",
    "\n",
    "a = BilinearForm(fes, eliminate_internal=condense)\n",
    "a += SymbolicBFI(grad(u)*grad(v))\n",
    "a += SymbolicBFI(alpha*order**2/h*jump_u*jump_v, element_boundary=True)\n",
    "a += SymbolicBFI(-grad(u)*n*jump_v - grad(v)*n*jump_u, element_boundary=True)\n",
    "\n",
    "b = CoefficientFunction( (20,1) )\n",
    "a += SymbolicBFI(-b * u * grad(v))\n",
    "uup = IfPos(b*n, u, uhat)\n",
    "a += SymbolicBFI(b*n*uup*jump_v, element_boundary=True)\n",
    "a.Assemble()\n",
    "\n",
    "f = LinearForm(fes)\n",
    "f += SymbolicLFI(1*v)\n",
    "f.Assemble()\n",
    "\n",
    "gfu = GridFunction(fes)\n",
    "\n",
    "print (\"A non-zero elements:\", a.mat.nze)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not condense:\n",
    "    inv = a.mat.Inverse(fes.FreeDofs(), \"umfpack\")\n",
    "    gfu.vec.data = inv * f.vec\n",
    "else:\n",
    "    f.vec.data += a.harmonic_extension_trans * f.vec \n",
    "    \n",
    "    inv = a.mat.Inverse(fes.FreeDofs(True), \"umfpack\")\n",
    "    gfu.vec.data = inv * f.vec\n",
    "    \n",
    "    gfu.vec.data += a.harmonic_extension * gfu.vec\n",
    "    gfu.vec.data += a.inner_solve * f.vec\n",
    "\n",
    "Draw (gfu.components[0], mesh, \"u-HDG\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
